Quasiperiodic patterns in boundary-modulated excitable waves 
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We investigate the impact of the domain shape on wave propagation in excitable media. Chan- 
nelled domains with sinusoidal boundaries are considered. Trains of fronts generated periodically at 
an extreme of the channel are found to adopt a quasiperiodic spatial configuration stroboscopically 
frozen in time. The phenomenon is studied in a model for the photo-sensitive Belousov-Zabotinsky 
reaction, but we give a theoretical derivation of the spatial return maps prescribing the height and 
position of the successive fronts that is valid for arbitrary excitable reaction-diffusion systems. 
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Excitable media display a very rich spatio-temporal 
behavior with regimes ranging from fairly well ordered 
structures of propagating waves (!]] to highly uncorrelated 
spatio-temporal chaos. The study of all these features as 
well as their mutual connections provides very useful in- 
sight to understand and eventually control phenomena of 
paramount applied importance such as the deadly arisal 
of fibrillation in cardiac tissues Q or the appearance of 
either ordered or turbulent patterns in extended chemi- 
cal reactors operating away from equilibrium conditions. 
In many of these applications, a crucial but frequently 
ignored ingredient is the presence of boundaries. For ex- 
ample, it has been shown that boundaries and obsta- 
cles in inhomogeneous media are important to either pin 
or repele spiral patterns 101; moving boundaries Q| and 
stripped configurations |j,|6| have also strong effects. 

Unfortunately, the current understanding of boundary 
effects in nonlinear partial differential equations is rather 
incomplete, and sometimes surprisingly nontrivial behav- 
ior lurk behind the apparent simplicity of some problems. 
A recent study (tJ] , for example, shows that relatively reg- 
ular boundary conditions such as Dirichlet's on the banks 
of a sausage-shaped channel can elicit several types of 
spatial complexity such as frozen quasiperiodicity and 
chaos even in very simple reaction diffusion equations. 
There, the axial coordinate along the channel acts as a 
"time" in the equations describing the time-independent 
spatial patterns and the undulated boundaries play the 
role of a periodic force inducing chaos in a dynamical 
system that is non-chaotic in the absence of driving. 

On the other hand, propagation of waves in excitable 
media has been studied in a variety of contexts (!]] . Due 
to their ubiquity in large two-dimensional systems, much 
of this work deals with spiral waves and focuses, in par- 
ticular, on the various aspects of the behavior around the 
cores of the spiral patterns. In contrast, the propagation 
of front trains has received much less attention. This 
may seem surprising since the same spirals can be seen 
far from their cores as a periodic train of two-dimensional 
traveling fronts. These trains, though, are easily charac- 



terized by a dispersion relation c = c(A), giving a rela- 
tion between the constant front train velocity and its 
uniform spacing A and their dynamics is very simple. 
However, much less trivial behavior appears even in one- 
dimensional systems if the excitable medium recovers the 
rest state not monotonically but via damped oscillations 
||. In this regime, propagating wave trains often relax 
to irregularly spaced configurations of fronts that can be 
seen as spatial chaos. 

The purpose of this Letter is to report a new kind 
of nontrivial spatial structure arising as a pure bound- 
ary effect in excitable media, namely, stroboscopically- 
frozen quasiperiodicity. We investigate the asymptotic 
propagation of excitable wave trains generated by lo- 
cal time-periodic stimulation at the extreme of a sinu- 
soidally undulated channel. We find that contrarily to 
intuitive expectations, the trains of fronts asymptotically 
accommodate in quasiperiodic spatial configurations, in- 
commensurated with the boundaries but periodic in time 
and synchronized with the stimuli. With the experiments 
on the photosensitive Belousov-Zhabotinsky reaction in 
mind , we demonstrate this phenomenon in the Oreg- 
onator model adapted to include the effect of light. Fi- 
nally, we present a more general semi-analytic theory of 
the formation of the quasi-periodic, and possibly chaotic, 
structures referred above. 

Photosensitive i?M(6&y)^" 2 -catalyzed Belousov- 



Zhabotinsky reactive media can be modelled 
following version of the Oregonator model: 



by the 



du 
~dt 



dv 

g~ t = (U ~ V) 



(fv + <f>) 



D u V 2 u 



(1) 



Here u (resp. v) describe HBr02 (resp. catalyst) con- 
centrations. D u is a diffusion coefficient and f,q,e and 
(j) are parameters related to the reaction kinetics. In our 
simulations we set / = 3, q = 0.002, e = 0.05, cf> = 0.002 
and D u = 1. 

We simulate this reaction in a spatial domain tailored 
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as an undulated channel along the longitudinal direction 
x. The transversal coordinate y is bounded by two si- 
nusoidal walls of common spatial frequency k = 2ir/\ p , 
amplitude d, minimum separation s, and phase mismatch 
77: 



Vo{x) = -[1 - cos(fcr)] 

y x (x) = -s - ~[1 + cos(fca; + 77)] 



(2) 



We will concentrate on symmetric sausage-shaped chan- 
nels (j) = 7r) of width w(x) = s + d — dcos(kx), as the 
one shown in Fig. [I]. On the sinusoidal boundaries we 
impose the Dirichlet condition u — 0.004, a value close 
to the fixed point of the local dynamics. 

In order to solve numerically Eq. ([j]) it is convenient 
to map the region limited by yo(x) and y\{x) and by 
x = 0, L to a rectangle: j/i = 1, yo = 0, and x = 0,Z, 
where y — (y — yo)/(yi — Uo) an d L is the length of the 
channel. Under this map, the diffusion term transforms 



V 2 u 



dt x u + F(x)dl~u + G(x)d z ~u - 

\ / yy \ / X y 



H(x)cHi, (3) 



F(x), G(x) and H(x), given in J7J, are periodic functions 
reflecting the undulations of the boundaries via modula- 
tions measured by the product kd. In the limit kd — > 
(straight channel), Eq.(0) becomes the standard Lapla- 



During the simulations we mainly varied the forcing pa- 
rameters A p and d, but also several wave train periods 
and channel widths were investigated. After a transient, 
the fields u and v converge to a configuration of prop- 
agating fronts that repeats itself periodically in time in 
synchrony with the wave generator at x = 0. In other 
words, the train becomes a stroboscopically frozen pat- 
tern. We denote by x n and a n the longitudinal position 
and maximum height of u at the channel axis, respec- 
tively, for the n-th front. 

As a comparison, in a straight channel (rf = 0) of finite 
width s the asymptotic configuration of the wave fronts 
is equally spaced by a length A and propagates with ve- 
locity c = X/T if the forcing period is T. This velocity 
increases with the channel width || starting from a crit- 
ical value s c of the latter, below which the fronts cannot 
propagate. In Fig. || we plot the train velocity and the 
maximum amplitude of the wave fronts as a function of 
the width of the straight channel. 
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FIG. 2. Wave train velocity and wave height (multiplied by 
a factor 10) in straight channels of different widths. There is 
a critical value, s c — 1.65, below which propagation becomes 
impossible. T = 5. 
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FIG. 1. In the upper panel, the white area limited by gray 
undulated boundaries is the excitable region where fronts orig- 
inated at the left end propagate. The transverse features are 
the fronts travelling to the right. Darker fronts have a larger 
value of u, as indicated in the colorbar. The lower panel dis- 
plays the envelope of the maximum amplitude of the fronts. 

Wave trains are generated stimulating the medium at 
the left end x — of the channel by pushing u above and 
below the excitability threshold periodically in time. The 
opposite end of the channel is set as a no-flux boundary. 



In modulated domains with d ^ a wide range of new 
spatial configurations incommensurated with the bound- 
aries emerge. Typically, both the spacing and the am- 
plitude of the fronts become spatially quasiperiodic. Ac- 
cording to the strength kd cx d/X p , of the spatial forcing 
we distinguish strong from weak modulations. Let us de- 
scribe the cases X p — 50 and X p = 1000, respectively, as 
an illustration. 

The results for strong modulation are shown in Fig. [jj. 
The amplitude of the boundary undulation increases 
from top to bottom. The quasiperiodic behavior of the 
pulse height becomes evident as d increases. The sec- 
ond column in Fig. ^| shows also the maximum a n of 
each front as a function of its position modulo A p . This 
plot provides information about the distribution of the 
front height maxima relative to the elementary unit of 
the channel. Notice that the fronts do not always reach 
their minimal height at the narrowest channel sections 
(x = mXp) as one would naively expect from the behav- 
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ior in straight channels depicted in Fig. ||. Moreover, the 
fronts can now propagate even when the channel is nar- 
rower (s = 1.1) in some places than the minimun width 
s c = 1.65 that allows propagation in straight channels. 
Last column in Fig. || displays the return maps of the 
(ri+ l)-th front position x n +i (relative to the unit channel 
cell) as a function of the position x n of the previous front. 
The shapes of the curves are analog to those circle maps 
describing the temporal dynamics of periodically forced 
sclf-oscillators. This analogy suggests that our system 
should exhibit the same richness of spatial behaviors as 
the circle map does in time-evolutions. 
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FIG. 3. Numerical results obtained from Eq. (jl|) for strong 
forcing. First column: Maximum height of each front within 
the train, as a function of position x. Second column: Same 
as before but with position x folded modulo \ v = 50. Third 
column: Return map of the front positions modulo \ v . Each 
row is for a different amplitude d. Parameters: A p = 50, 
T — 5, s — 1.1. Note that channel length L = 15000 is much 
larger that in Fig. jjj. 

The weak forcing case is illustrated in Fig. ^. As in the 
case of circle maps for very weak forcing, the front po- 
sitions return map shows a very small deviation from 
linearity with the given parameter values. This approx- 
imate linearity implies that the front train wavelength 
is nearly constant and the influence of the channel walls 
negligible. This influence is, however, more important 
on the front heights. Minima of front height are situated 
at the narrowest channel sections, in concordance with 
Fig. |^, while the maxima saturate for d large enough. 
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FIG. 4. Weak forcing behavior at A p = 1000 and s — 2. 
Above: Return map of the front positions modulo X p . Be- 
low: Maximum height of each front as a function of the front 
position modulo X p . Solid lines are from Eqs. (|) and @. 

Let us now derive a semi- analytical expression for the 
return maps of successive wave fronts positions and max- 
ima heights. In view of the results for the weak forcing 
case we assume that the front velocity in the undulated 
channel at a position where the local width is w adapts 
quasi-adiabatically to the velocity c(w) (Fig. ||) corre- 
sponding to a uniform channel of the same width. Thus, 
the velocity of the n-th front is 



x n (t) = c(w(x n )) 



In our channel w(x) = wq — wi cos(fcx) with wq = s + d 
and w\ = d. In order to proceed analytically an approx- 
imation for c(w) should be introduced. For d small the 
width variation is also small and c(w) can be replaced by 
a linear fit a + bw of an apropriate range of data in Fig. |^. 
Hence, c(w(x n )) « cq — c\ cos(kx n ) where cq — a + bwo 
and ci = bw\. Eq. (|J) can now be integrated during one 
period T of the front generator, to obtain: 



k^/cj 



[arct an /(£„)] 



x„ +I (0) 

x n (0) 



= T 



(5) 



with 



/On) 



CO 



c - ci 



Ci / kx n 
— tan — 
V 2 



zrii (6) 



Here we have used the observed time periodicity of the 
wave train to write x n (T) = a; n +i(0), a crucial step 
to convert the time-differential equation (^) into a map 
for space positions. Defining ip = arctanz and A = 
0.5 kTyJcq — c\ we have A = ip n +i — ip n , and the return 
map for the variable z is 

Zn+l = g{Zn) 



tan A 



1 — tan A • z n 
In terms of the front position x we finally have: 



(7) 
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X n+1 = f- 1 (9(f(Xn))) (8) 

For the maximum height of the wave fronts, 
the same adiabaticity assumption leads to a n — 
h (wq ~ wi cos(kx n )), with h(w) being the maximum 
height of the fronts in a straight channel of width w. We 
can go one step further towards qualitatively describe 
the observed positional mismatch between the minimal- 
height fronts and the narrowest channel sections by con- 
sidering a short adaptation time r a of the front charac- 
teristics to the local width: 

a n {t) — — (h(w(x n )) - an) (9) 

Ta 

As above, by linearly fitting the data from Fig. || in the 
range (s, s + 2d) we have h(w) ss a' + b'w for small d. 
Then, h(w(x)) « ho — hi cos kx, with ho = a' + b'wo, and 
hi = b'wi. Integrating Eq. (^) for small d and kcoT a <C 1 
so that we can set x n {t) — a;„(0) + cot + 0(d), we get a 
relationship linking the wavefront heights and positions, 

a n +i = h , hl n = sin (k(x n + v T) + 6) . (10) 

Here 9 — arccos(l/yl -I- (r a vok) 2 ) describes the dis- 
placement of the minimal heights from the narrowest sec- 
tions. 



our boundary-induced patterns in a given limit. Finally, 
while the agreement between the theory and the numerics 
is bound to worsen as forcing increases, the theory still 
describes well the pattern features in the strong forcing 
regime. For instance, Fig. |B| shows how the maxima and 
minima of a n (x n ) shift as d is increased, in good quali- 
tative agreement with the corresponding case of Fig. ||. 

In summary, we have shown that boundary conditions 
in domains with the form of undulated channels may 
induce nontrivial longitudinal spatial configurations of 
propagating trains of excitation fronts generated by a 
local time-periodic stimulation in simple excitable me- 
dia. In particular, stroboscopically frozen quasi-periodic 
arrays of fronts were found. These structures were de- 
scribed in terms of spatial return maps that are very 
similar to the circle maps whose iteration describe the 
temporal dynamics of forced oscillators. This similar- 
ity allows one to speculate about the existence of even 
more complex configurations representing the spatial re- 
alizations of the chaotic regimes of these maps. The phe- 
nomenon reported here should be experimentally observ- 
able in the photo-sensitive Belusov-Zhabotinsky reaction 
with proper lighting conditions at the boundaries. Work 
along this line is currently in progress. 
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0.6 p 
0.5 
0.4 
0.3 

0.2° 
0.1"° 



d=l 
d=2 



10 



20 



30 



40 



50 



x mod ( X ) 
n p 



FIG. 5. Change in the shape of a„(x„) from Eq. ( [lo| ) by 
increasing d beyond the weak forcing limit. The qualitative 
trend agrees with Fig. |j| 



Since the derivation of Eqs. (g) and ( |10| ) is formally valid 
only in the weak forcing limit we first contrast the the- 
ory against the numerical data in Fig. || for d = 0.5, to 
confirm the good agreement |I0[ . More detailed numer- 
ical explorations reasure us that both adiabaticity and 
small d approximations are justified and that the small 
deviations in Fig. |3| are only due to the linear approxi- 
mation. Moreover, a systematic d-expansion in Eq. (j^) 
would lead precisely to a circle map supporting the ob- 
servation that this model is relevant to the description of 
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